Questions abordées dans ce cours

1. Les données, leurs représentations

• data BIC issue de la base TCGA

2. Comment comparer deux individus

• notion de distance

3. Comment découvrir des “clusters” dans les données ?

• classification hiérarchique
• kmeans

4. Comment déterminer le nombre de groupe optimal ?
5. Comment comparer deux classifications ?

Les données : TCGA (1)

Les données sont issues de la base Recount2. Nous avons sélectionné l’étude TCGA : The Cancer Genome Atlas, regroupant des données RNA-seq pour plus de 12.000 patients souffrant de différents types de cancer. Nous nous intéressons ici uniquement aux données Breast Invasive Cancer (BIC) concernant le cancer du sein.

Les données ont été préparées pour vous, selon la procédure détaillée au cours sur l’analyse différentielle de données RNA-seq.

1. Filtrage des gènes à variance nulle et de ceux contenant trop de zéros.

2. Normalisation (méthode robuste aux outliers)

3. Analyse différentielle multi-groupes (en utilisant le package Bioconductor edgeR).

4. Correction des P-valeurs pour tenir compte des tests multiples (nous avons testé ici ~20.000 gènes). Nous estimons le False Discovery Rate (FDR) selon la méthode de Benjamini-Hochberg (fonction R p.adjust(all.pvalues, method="fdr")).

5. Sélection de gènes différentiellement exprimés sur base d’un seuil \(\alpha = 0.05\) appliqué au FDR.

TCGA (2)

         637E.4A4 11F.41D3 7FB8.471 EC6.4FB2
00003.14 16.52775 18.84599 18.67677 17.75458
00419.12 17.50349 17.01692 18.15315 19.34170
00457.13 17.96040 17.29379 17.51914 18.34039
00460.16 16.77151 16.30180 16.85386 17.75886
00938.12 15.13011 13.57196 15.04359 15.11192
00971.15 17.21899 20.05815 17.48420 18.33249

Pour des raisons historiques, en analyse transcriptomique les données sont toujours fournies avec

• 1 ligne = 1 gène
• 1 échantillon biologique = 1 colonne

Cette convention a été établie en 1997, lors des toutes premières publications sur le transcriptome de la levure. Dans ces études, l’objet d’intérêt (l’“individu”) était le gène, et les variables étaient ses mesures d’expression dans les différentes conditions testées.

Pour l’analyse de tissus cancéreux, on considère au contraire que l’“objet” d’intérêt est l’échantillon prélevé sur le patient, et les variables sont les mesures d’expression des différents gènes chez un patient.

TCGA (3)

Classiquement, en analyse de données, les individus sont les lignes du tableau de données, les colonnes sont les variables.

Ce qui implique de faire attention, et éventuellement de travailler sur la matrice transposée (fonction t() en R) pour utiliser correctement les fonctions classiques.

t(BIC.expr[1:6, 1:4])

         00003.14 00419.12 00457.13 00460.16 00938.12 00971.15
637E.4A4 16.52775 17.50349 17.96040 16.77151 15.13011 17.21899
11F.41D3 18.84599 17.01692 17.29379 16.30180 13.57196 20.05815
7FB8.471 18.67677 18.15315 17.51914 16.85386 15.04359 17.48420
EC6.4FB2 17.75458 19.34170 18.34039 17.75886 15.11192 18.33249

• 1 ligne = 1 gène = 1 individu = 1 vecteur

• 1 colonne = 1 feature = 1 vecteur

• l’ensemble des données = 1 data.frame

[1]  819 1000

Représentons ces données (1)

• extrait des données
  • chaque individu est représenté par un vecteur de mesures

         00003.14 00419.12 00457.13 00460.16 00938.12
637E.4A4 16.52775 17.50349 17.96040 16.77151 15.13011
11F.41D3 18.84599 17.01692 17.29379 16.30180 13.57196
7FB8.471 18.67677 18.15315 17.51914 16.85386 15.04359
EC6.4FB2 17.75458 19.34170 18.34039 17.75886 15.11192

• Comment représenter / visualiser ces données ?
• Dans quel espace de réprésentation ?

Représentons ces données : un individu à la fois (2)

Représentons ces données : deux individu à la fois (3)

Représentons ces données : une variable à la fois (4)

Représentons ces données : deux variables à la fois (5)

Représentons ces données : toutes les variables (6)

en tenant compte de l’ensemble des individus/ lignes et variables / colonnes = un nuage de points dans un espace à 1000 dimensions

• chaque point est représenté par un vecteur dans \(\mathbb{R}^{1000}\)
• le nuage de points est représenté par une matrice à n (= 819) et p (= 1000 dimensions)
  • n = nombre de lignes = nombre d’individus = taille de l’échantillon
  • p = nombre de colonnes = nombre de variables décrivant l’échantillon

= PAS de représentation possible (pour l’instant)

On a cependant

• ACP
• heatmap

ACP

heatmap

Clustering et classification : définition

Nous utiliserons les termes anglais

en français :

• clustering = classification
• classification = discrimination

Clustering

On a pas d’information supplémentaire sur nos données, juste le data.frame contenant

• variables quantitatives = vecteur de réels

Clustering : on cherche à mettre en évidence des groupes (/ des clusters) dans les données

• un groupe =
  • des individus qui se ressemblent et
  • qui sont différents des autres groupes
• le clustering appartient aux méthodes dites non supervisées, ou descriptives

Classification

On a une information** supplémentaire sur nos données : on connaît le partitionnement de notre jeu de données

• variables quantitatives = vecteur de réels
• ET
• variable qualitative = groupe (cluster)

Classification : on cherche un algorithme / un modèle permettant de prédire la classe, le groupe de tout individu dont on connait les caractéristiques

• la classification appartient aux méthodes dites supervisées, ou prédictives

Clustering

- y a t’il des groupes ? si oui, combien ?

• Méthode agglomérative ou hierarchical clustering
• Moyennes mobiles ou K-means : séparation optimale des groupes connaissant le nombre de groupes

Comment comparer des vecteurs-individus ?

Distances

Définition d’une distance : fonction positive de deux variables

1. \(d(x,y) \ge 0\)
2. \(d(x,y) = d(y,x)\)
3. \(d(x,y) = 0 \Longleftrightarrow x = y\)
4. Inégalité triangulaire : \(d(x,z) \le\) d(x,y)+d(y,z)

Si 1,2,3 seulement: dissimilarité

Distance euclidienne et distance de corrélation

Avec R (1) : distance entre deux individus

• on utilise la fonction dist() avec l’option method = "euclidean", "manhattan", ...

distance euclidienne : dist(mat.xy) = 2.38

distance de manhattan = dist(mat.xy, method = "manhattan") 4.38

• on utilise la fonction 1 - cor() avec l’option method = "pearson", "spearman", ...

distance de corrélation = 1-cor(t(mat.xy) 0.75

Avec R (2) : distance entre individus d’un nuage de points

• distance euclidienne, de 5 individus choisis au hasard

         2D80.49F D91.450F 0DD.4F9A 573E.42E
D91.450F       35                           
0DD.4F9A       33       33                  
573E.42E       53       54       53         
CE60.480       31       29       23       52

• distance de corrélation : \(d = 1-r\)

         2D80.49F D91.450F 0DD.4F9A 573E.42E
D91.450F    0.149                           
0DD.4F9A    0.143    0.140                  
573E.42E    0.287    0.293    0.294         
CE60.480    0.118    0.111    0.076    0.279

Avec R (3) : distance entre variables décrivant le nuage de points

         00003.14 00419.12 00457.13 00460.16 00938.12
00419.12    1.270                                    
00457.13    0.615    1.739                           
00460.16    0.651    1.402    0.087                  
00938.12    1.771    0.463    1.381    1.162         
00971.15    0.458    1.875    0.112    0.335    1.731

Distances entre groupes (1)

Distances entre groupes (2)

Les données

Revenons à nos données

Visualisation des données (1)

On peut ensuite essayer de visualiser les données

• par un plot (! ne pas faire si “grosses” données)

Visualisation des données (2)

• par un boxplot (! ne pas faire si “grosses” données)

boxplot(BIC_expr[, sample(1:ncol(BIC_expr), 30)], las = 2)

Préparation des données (1) : variables de variance nulle

BIC.var <- apply(BIC_expr, 2, var)
sum(apply(BIC_expr, 2, var) == 0)

[1] 0

Préparation des données (2) : Mise à l’échelle

Afin de pouvoir considérer que toutes les variables sont à la même échelle, il est parfois nécessaire de standardiser les données.

• soit

  • en centrant (ramener la moyenne / médiane de chaque variable à \(0\))

BIC.expr.centre <- scale(BIC_expr, center = TRUE, scale = FALSE)

• soit

  • en centrant (ramener la moyenne de chaque variable \(0\))
  • et mettant à l’échelle (ramener la variance de chaque variable à \(1\))

BIC.expr.scaled <- scale(BIC_expr, center = TRUE, scale = TRUE)

• soit en effectuant une transformation des variables, par exemple transformation logarithmique

On peut visuellement regarder l’effet de la standardisation

• par des boîtes à moustaches (boxplot)

Centrage sur la moyenne ou la médiane

Mise à l’échelle écart-type ou intervalle interquartile

Standardisation : centrage et mise à l’échelle

La classification hiérarchique : principe

classification hiérarchique : mettre en évidence des liens hiérachiques entre les individus

• classification hiérarchique ascendante : partir des individus pour arriver à des classes / cluster
• classification hiérarchique descendante : partir d’un groupe qu’on subdivise en sous-groupes /clusters jusqu’à arriver à des individus.

Notion importante, cf distances

• ressemblance entre individus = distance

  • euclidienne
  • corrélation

• ressemblance entre groupes d’invidus = critère d’aggrégation

  • lien simple
  • lien complet
  • lien moyen
  • critère de Ward

L’algorithme : étape 1

• départ : n individus = n clusters distincts

• calcul des distances entre tous les individus

  • choix de la métrique à utiliser en fonction du type de données

• regroupement des 2 individus les plus proches => (n-1) clusters

Au départ

Identification des individus les plus proches

Construction du dendrogramme

Etape j :

• calcul des dissemblances entre chaque groupe obtenu à l’étape \((j-1)\)

• regroupement des deux groupes les plus proches => \((n-j)\) clusters

Calcul des nouveaux représentants ‘BE’ et ‘CD’

Calcul des distances de l’individu restant ‘A’ aux points moyens

A est plus proche de …

dendrogramme

pour finir

• à l’étape \((n-1)\), tous les individus sont regroupés dans un même cluster

dendrogramme final

Sur nos données (1) : deux métriques différentes

Sur nos données (2) : deux critères d’aggrégation différents

En conclusion

• Faire attention au données

  • données manquantes
  • données invariantes
  • données normalisées

• Choisir

  • la distance (entre individus) et
  • le critère d’aggrégation (entre cluster)

adaptés à nos données

Visualisation à l’aide de heatmap : données brutes

pheatmap::pheatmap(BIC_expr, clustering.method = "ward.D2")

Visualisation à l’aide de heatmap : données mise à l’échelle

pheatmap::pheatmap(BIC.expr.scaled, clustering.method = "ward.D2")

Les k-means

Les individus dans le plan

=> faire apparaitres des classes / des clusters

L’algorithme

étape 1 :

• \(k\) centres provisoires tirés au hasard
• \(k\) clusters créés à partir des centres en regroupant les individus les plus proches de chaque centre
• obtention de la partition \(P_0\)

Choix des centres provisoires

Calcul des distances aux centres provisoires

Affectation à un cluster

Calcul des nouveaux centres de classes

Etape j :

• construction des centres de gravité des k clusters construits à l’étape \((j-1)\)

• \(k\) nouveaux clusters créés à partir des nouveaux centres suivant la même règle qu’à l’étape \(0\)

• obtention de la partition \(P_j\)

Fin :

• l’algorithme converge vers une partition stable

Arrêt :

• lorsque la partition reste la même, ou lorsque la variance intra-cluster ne décroit plus, ou lorsque le nombre maximal d’itérations est atteint.

Un premier k-means en 5 groupes

iris.scale.kmeans5 <- kmeans(BIC.expr.scaled, center=5)
iris.scale.kmeans5$cluster[1:20]

637E.4A4 11F.41D3 7FB8.471 EC6.4FB2 713.4383 9AA.4949 9951.4D4 AF7.4CBC EAD9.4BC 0319.4B5 A84.4B0B FDC1.414 008.44AF D0B.49F1 0930.403 3507.496 613.45DA C250.4F5 8AE.4245 792E.4DA 
       1        5        3        3        5        5        1        1        1        3        1        5        5        2        1        1        1        2        1        1 

Comment déterminer le nombre de clusters ? (1)

Ces méthodes non supervisées, sont sans a priori sur la structure, le nombre de groupe, des données.

rappel : un cluster est composé

• d’individus qui se ressemblent
• d’individus très différents des individus de ceux des autres clusters

Comment déterminer le nombre de clusters ? (2)

• si les individus d’un même cluster sont proches

  • homogénéité maximale à l’intérieur de chaque cluster => variance intra faible

• si les individus de 2 clusters différents sont éloignés => variance inter forte

  • hétérogénéité maximale entre chaque cluster

Comment déterminer le nombre de clusters ? avec la classification hiérarchique

La coupure de l’arbre à un niveau donné construit une partition. la coupure doit se faire :

• après les agrégations correspondant à des valeurs peu élevées de l’indice

• avant les agrégations correspondant à des niveaux élevés de l’indice, qui dissocient les groupes bien distincts dans la population.

Comment déterminer le nombre de clusters ? avec les kmeans

Comparaison des résultats des deux clustering

• par une table

Pros et cons des différents algorithmes

Comparaison de clustering: Rand Index

Mesure de similarité entre deux clustering

à partir du nombre de fois que les clustering sont d’accord

\[R=\frac{m+s}{t}\]

• \(m\) = nombre de paires dans la même classe dans les deux classifications
• \(s\) = nombre de paires séparées dans les deux classifications
• \(t\) = nombre total de paires

## Compute Rand index
(RI <- aricode::RI(cluster.hclust4, cluster.kmeans3))

[1] 0.705682

Comparaison de clustering: Adjusted Rand Index

\[ \text{ARI} = \frac{\text{RI}-\text{E(RI)}}{\text{Max RI} - \text{E(RI)}}\]

• \(\text{ARI}\) = adjusted Rand Index = RI normalisé
• \(E(RI)\) = expected RI, espérance aléatoire (en assignant les groupes au hasard)
• Prend en compte la taille des classes
• \(\text{ARI}= 1\) pour classification identique
• \(\text{ARI} \simeq 0\) pour classification aléatoire (peut être <0)
• Adapté même si les nombres de classes diffèrent entre les deux classifications
• Adapé à des tailles de classes différentes

Comparaison des résultats des deux classifications

• rand index et adjusted rand index

## Compute adjusted Rand index
(ARI <- aricode::ARI(cluster.hclust4, cluster.kmeans3))

[1] 0.3734816

Supplementary materials

POUR ALLER PLUS LOIN

Distances utilisées dans R (1)

• distance euclidienne ou distance \(L_2\): \(d(x,y)=\sqrt{\sum_i (x_i-y_i)^2}\)

• distance de manahattan ou distance \(L_1\): \(d(x,y)=\sum_i |x_i-y_i|\)

• distance du maximum ou L-infinis, \(L_\infty\): \(d(x,y)=\max_i |x_i-y_i|\)

Distances utilisées dans R (2)

• distance de Minkowski \(l_p\): \[d(x,y)=\sqrt[p]{\sum_i (|x_i-y_i|^p}\]

• distance de Canberra (x et y valeurs positives): \[d(x,y)=\sum_i \frac{x_i-y_i}{x_i+y_i}\]

• distance binaire ou distance de Jaccard ou Tanimoto: proportion de propriétés communes

Note : lors du TP, sur les données d’expression RNA-seq, nous utiliserons le coefficient de corrélation de Spearman et la distance dérivée, \(d_c = 1-r\)

Autres distances non géométriques (pour information)

Utilisées en bio-informatique:

• Distance de Hamming: nombre de remplacements de caractères (substitutions)

• Distance de Levenshtein: nombre de substitutions, insertions, deletions entre deux chaînes de caractères

\[d("BONJOUR", "BONSOIR")=2\]

• Distance d’alignements: distances de Levenshtein avec poids (par ex. matrices BLOSSUM)

• Distances d’arbre (Neighbor Joining)

• Distances ultra-métriques (phylogénie UPGMA)

Distances plus classiques en génomique

Il existe d’autres mesures de distances, plus ou moins adaptées à chaque problématique :

• Jaccard (comparaison d’ensembles): \(J_D = \frac{A \cap B}{A \cup B}\)

• Distance du \(\chi^2\) (comparaison de tableau d’effectifs)

Ne sont pas des distances, mais indices de dissimilarité :

• Bray-Curtis (en écologie, comparaison d’abondance d’espèces)

• Jensen-Shannon (comparaison de distributions) # Distance avec R : indice de Jaccard

• ou pour des distances particulières, par exemple l’indice de Jaccard :

          v.a       v.b
v.b 0.3333333          
v.c 0.0000000 0.3333333

… par une projection sur une ACP

par(mfrow = c(1,2))
biplot(prcomp(BIC.expr), las = 1, cex = 0.7,
       main = "Données non normalisées")
biplot(prcomp(BIC.expr, scale = TRUE), las = 1, cex = 0.7,
       main = "Données normalisées")

Géométrie et distances (1)

On considère les données comme des points de \(\mathbb{R}^n\)

\(\mathbb{R}^n\) : espace Euclidien à \(n\) dimensions, où

• chaque dimension représente une des variables observées;
• un individu est décrit comme un vecteur à \(n\) valeurs, qui correspond à un point dans cet espace.

Géométrie et distances (2)

On considère les données comme des points de \(R^n\) (*)

• géométrie donnée par distances
• distances = dissimilarités imposées par le problème
• dissimilarités \(\longrightarrow\) permettent la visualisation de l’ensemble des points

Géométrie et distances (3)

Sur la base d’une distance

• Clustering :

## Plot distances between 3 points in a 2D Euclidian space
plot(x = 0, y = 0, type = "n", xlim = c(0, 5), ylim = c(0, 5), 
    xlab = "", ylab = "", las = 1, 
     main = "3 individuals in a 2-D space\nDot plot representation",
     panel.first = grid())
points(x = 1, y = 1, col = "blue", pch = 19)
text(x = 1, y = 1, col = "blue", label = "A", pos = 2)
points(x = 2, y = 0, col = "green", pch = 19)
text(x = 2, y = 0, col = "green", label = "B", pos = 4)
points(x = 4, y = 4, col = "red", pch = 19)
text(x = 4, y = 4, col = "red", label = "C", pos = 4)

Distance euclidienne

• distance euclidienne ou distance \(L_2\): \(d(x,y)=\sqrt{\sum_i (x_i-y_i)^2}\)

Distances entre groupes (2)

• Single linkage : élements les plus proches des 2 groupes

\[D(C_1,C_2) = \min_{i \in C_1, j \in C_2} D(x_i, x_j)\]

• Complete linkage : éléments les plus éloignés des 2 groupes

\[D(C_1,C_2) = \max_{i \in C_1, j \in C_2} D(x_i, x_j)\]

• Average linkage : distance moyenne

\[D(C_1,C_2) = \frac{1}{N_1 N_2} \sum_{i \in C_1, j \in C_2} D(x_i, x_j)\]

• Ward

\(d^2(C_i,C_j) = I_{intra}(C_i \cup C_j)-I_{intra}(C_i)-I_{intra}(C_j)\)

\(D(C_1,C_2) = \sqrt{\frac{N_1N_2}{N_1 + N_2}} \| m_1 -m_2 \|\)

La matrice de distance euclidienne

La matrice de distance de corrélation

Les heatmap - échelle de couleur standardisée par colonne

pheatmap::pheatmap(BIC_expr, scale = "column", clustering.method = "ward.D2")

Les heatmap - échelle de couleur standardisée par ligne

pheatmap::pheatmap(BIC_expr, scale = "row", clustering.method = "ward.D2")

Avec R (1) : distance entre deux individus

• on utilise la fonction dist() avec l’option method = "euclidean", "manhattan", ...

distance euclidienne : 2.91

distance de manhattan = 5.35

• on utilise la fonction 1 - cor() avec l’option method = "pearson", "spearman", ...

distance de corrélation = 1.95

Je ne fais pas attention à ce que je fais …

… c’est à dire aux options des fonctions dist() et hclust()

par(mfrow = c(2, 1))
plot(BIC.hclust, hang = -1, cex = 0.5, main = "Données brutes")
plot(BIC.scale.hclust, hang = -1, cex = 0.5, main = "Normalisées")

Supplément : analyse de données d’expression 2019

• TP clustering : [html] [pdf] [Rmd]

• Première partie : chargement des données

Contact: anne.badel@univ-paris-diderot.fr